
#ifndef GH_MATH_H_INCLUDED
#define GH_MATH_H_INCLUDED

#include <stdint.h>
#include <math.h>
#include <xmmintrin.h>
#include <emmintrin.h>

#ifndef GH_INLINE
    #if defined(_MSC_VER)
        #define GH_INLINE __forceinline
    #else
        #define GH_INLINE inline
    #endif
#endif

#if defined(_MSC_VER)
    #pragma warning(disable: 4201) /* nameless struct/union */
#endif

#ifndef GH_ALIGN
    #if defined(_MSC_VER)
        #define GH_ALIGN(x) __declspec(align(x))
    #else
        #define GH_ALIGN(x) __attribute__((aligned(x)))
    #endif
#endif

/* Math Constants */
#define GH_PI 3.141592654f
#define GH_2PI 6.283185307f
#define GH_1DIVPI 0.318309886f
#define GH_1DIV2PI 0.159154943f
#define GH_PIDIV2 1.570796327f
#define GH_PIDIV4 0.785398163f

#define GH_SHUFFLE(x, y, z, w) _MM_SHUFFLE(w, z, y, x)
#define GH_X 0
#define GH_Y 1
#define GH_Z 2
#define GH_W 3
#define ghVecShuffle(v0, v1, mask) _mm_shuffle_ps((v0), (v1), (mask))

/* GHvec */
typedef __m128 GHvec;

/* GHvecF32 */
typedef GH_ALIGN(16) struct {
    union {
        float f[4];
        GHvec v;
    };
} GHvecF32;

/* GHvecS32 */
typedef GH_ALIGN(16) struct {
    union {
        int32_t s[4];
        GHvec v;
    };
} GHvecS32;

/* GHvecU32 */
typedef GH_ALIGN(16) struct {
    union {
        uint32_t u[4];
        GHvec v;
    };
} GHvecU32;

/* GHmat */
typedef GH_ALIGN(16) struct {
    union {
        GHvec r[4];
        struct {
            float _00, _01, _02, _03;
            float _10, _11, _12, _13;
            float _20, _21, _22, _23;
            float _30, _31, _32, _33;
        };
        float m[4][4];
    };
} GHmat;

/* GHfloat2 */
typedef struct {
    float x;
    float y;
} GHfloat2;

/* GHfloat3 */
typedef struct {
    float x;
    float y;
    float z;
} GHfloat3;

/* GHfloat4 */
typedef struct {
    float x;
    float y;
    float z;
    float w;
} GHfloat4;

/* GHfloat3x3 */
typedef struct {
    union {
        struct {
            float _00, _01, _02;
            float _10, _11, _12;
            float _20, _21, _22;
        };
        float m[3][3];
    };
} GHfloat3x3;

/* GHfloat4x3 */
typedef struct {
    union {
        struct {
            float _00, _01, _02;
            float _10, _11, _12;
            float _20, _21, _22;
            float _30, _31, _32;
        };
        float m[4][3];
    };
} GHfloat4x3;

/* GHfloat4x4 */
typedef struct {
    union {
        struct {
            float _00, _01, _02, _03;
            float _10, _11, _12, _13;
            float _20, _21, _22, _23;
            float _30, _31, _32, _33;
        };
        float m[4][4];
    };
} GHfloat4x4;

/* Aligned Versions */
typedef GH_ALIGN(16) GHfloat2 GHfloat2A;
typedef GH_ALIGN(16) GHfloat3 GHfloat3A;
typedef GH_ALIGN(16) GHfloat4 GHfloat4A;
typedef GH_ALIGN(16) GHfloat4x3 GHfloat4x3A;
typedef GH_ALIGN(16) GHfloat4x4 GHfloat4x4A;

/* FloatN Set Operations */
static GHfloat2 ghFloat2Set(float x, float y);
static GHfloat3 ghFloat3Set(float x, float y, float z);
static GHfloat4 ghFloat4Set(float x, float y, float z, float w);

static GHfloat2 ghFloat2Setv(const float *src);
static GHfloat3 ghFloat3Setv(const float *src);
static GHfloat4 ghFloat4Setv(const float *src);

static GHfloat2A ghFloat2ASet(float x, float y);
static GHfloat3A ghFloat3ASet(float x, float y, float z);
static GHfloat4A ghFloat4ASet(float x, float y, float z, float w);

static GHfloat2A ghFloat2ASetv(const float *src);
static GHfloat3A ghFloat3ASetv(const float *src);
static GHfloat4A ghFloat4ASetv(const float *src);

/* FloatMxN Set Operations */
static GHfloat3x3 ghFloat3x3Set(float, float, float, float, float, float, float, float, float);
static GHfloat4x3 ghFloat4x3Set(float, float, float, float, float, float, float, float, float, float, float, float);
static GHfloat4x4 ghFloat4x4Set(float, float, float, float, float, float, float, float, float, float, float, float, float, float, float, float);

static GHfloat3x3 ghFloat3x3Setv(const float *src);
static GHfloat4x3 ghFloat4x3Setv(const float *src);
static GHfloat4x4 ghFloat4x4Setv(const float *src);

static GHfloat4x3A ghFloat4x3ASet(float, float, float, float, float, float, float, float, float, float, float, float);
static GHfloat4x4A ghFloat4x4ASet(float, float, float, float, float, float, float, float, float, float, float, float, float, float, float, float);

static GHfloat4x3A ghFloat4x3ASetv(const float *src);
static GHfloat4x4A ghFloat4x4ASetv(const float *src);

/* Vector Load Operations */
static GHvec ghVecLoadFloat(const float *src);
static GHvec ghVecLoadFloat2(const GHfloat2 *src);
static GHvec ghVecLoadFloat3(const GHfloat3 *src);
static GHvec ghVecLoadFloat4(const GHfloat4 *src);
static GHvec ghVecLoadFloat2A(const GHfloat2A *src);
static GHvec ghVecLoadFloat3A(const GHfloat3A *src);
static GHvec ghVecLoadFloat4A(const GHfloat4A *src);

/* Matrix Load Operations */
static GHmat ghMatLoadFloat3x3(const GHfloat3x3 *src);
static GHmat ghMatLoadFloat4x3(const GHfloat4x3 *src);
static GHmat ghMatLoadFloat4x4(const GHfloat4x4 *src);
static GHmat ghMatLoadFloat4x3A(const GHfloat4x3A *src);
static GHmat ghMatLoadFloat4x4A(const GHfloat4x4A *src);

/* Vector Store Operations */
static void ghVecStoreFloat(float *dst, GHvec v);
static void ghVecStoreFloat2(GHfloat2 *dst, GHvec v);
static void ghVecStoreFloat3(GHfloat3 *dst, GHvec v);
static void ghVecStoreFloat4(GHfloat4 *dst, GHvec v);
static void ghVecStoreFloat2A(GHfloat2A *dst, GHvec v);
static void ghVecStoreFloat3A(GHfloat3A *dst, GHvec v);
static void ghVecStoreFloat4A(GHfloat4A *dst, GHvec v);

/* Matrix Store Operations */
static void ghMatStoreFloat3x3(GHfloat3x3 *dst, const GHmat *m);
static void ghMatStoreFloat4x3(GHfloat4x3 *dst, const GHmat *m);
static void ghMatStoreFloat4x4(GHfloat4x4 *dst, const GHmat *m);
static void ghMatStoreFloat4x3A(GHfloat4x3A *dst, const GHmat *m);
static void ghMatStoreFloat4x4A(GHfloat4x4A *dst, const GHmat *m);

/* Vector Set and Get Operations */
static GHvec ghVecZero(void);
static GHvec ghVecOne(void);
static GHvec ghVecEpsilon(void);
static GHvec ghVecSet(float x, float y, float z, float w);
static GHvec ghVecReplicate(float s);
static GHvec ghVecReplicatev(const float *src);
static float ghVecGetX(GHvec v);
static float ghVecGetY(GHvec v);
static float ghVecGetZ(GHvec v);
static float ghVecGetW(GHvec v);

/* Integer Vector Set Operations */
static GHvec ghVecITrue(void);
static GHvec ghVecIFalse(void);
static GHvec ghVecISet(uint32_t x, uint32_t y, uint32_t z, uint32_t w);
static GHvec ghVecIReplicate(uint32_t s);
static GHvec ghVecIReplicatev(const uint32_t *src);

/* Vector Comparison Operations */
static GHvec ghVecEqual(GHvec v0, GHvec v1);
static GHvec ghVecNearEqual(GHvec v0, GHvec v1, GHvec epsilon);
static GHvec ghVecNotEqual(GHvec v0, GHvec v1);
static GHvec ghVecLess(GHvec v0, GHvec v1);
static GHvec ghVecLessEqual(GHvec v0, GHvec v1);
static GHvec ghVecGreater(GHvec v0, GHvec v1);
static GHvec ghVecGreaterEqual(GHvec v0, GHvec v1);

/* VectorN Comparison Operations */
static int ghVec2Equal(GHvec v0, GHvec v1);
static int ghVec3Equal(GHvec v0, GHvec v1);
static int ghVec4Equal(GHvec v0, GHvec v1);
static int ghVec2NearEqual(GHvec v0, GHvec v1, GHvec epsilon);
static int ghVec3NearEqual(GHvec v0, GHvec v1, GHvec epsilon);
static int ghVec4NearEqual(GHvec v0, GHvec v1, GHvec epsilon);
static int ghVec2NotEqual(GHvec v0, GHvec v1);
static int ghVec3NotEqual(GHvec v0, GHvec v1);
static int ghVec4NotEqual(GHvec v0, GHvec v1);
static int ghVec2Less(GHvec v0, GHvec v1);
static int ghVec3Less(GHvec v0, GHvec v1);
static int ghVec4Less(GHvec v0, GHvec v1);
static int ghVec2LessEqual(GHvec v0, GHvec v1);
static int ghVec3LessEqual(GHvec v0, GHvec v1);
static int ghVec4LessEqual(GHvec v0, GHvec v1);
static int ghVec2Greater(GHvec v0, GHvec v1);
static int ghVec3Greater(GHvec v0, GHvec v1);
static int ghVec4Greater(GHvec v0, GHvec v1);
static int ghVec2GreaterEqual(GHvec v0, GHvec v1);
static int ghVec3GreaterEqual(GHvec v0, GHvec v1);
static int ghVec4GreaterEqual(GHvec v0, GHvec v1);

/* Integer Vector Comparison Operations */
static GHvec ghVecIEqual(GHvec v0, GHvec v1);
static GHvec ghVecINotEqual(GHvec v0, GHvec v1);

/* Integer VectorN Comparison Operations */
static int ghVec2IEqual(GHvec v0, GHvec v1);
static int ghVec3IEqual(GHvec v0, GHvec v1);
static int ghVec4IEqual(GHvec v0, GHvec v1);
static int ghVec2INotEqual(GHvec v0, GHvec v1);
static int ghVec3INotEqual(GHvec v0, GHvec v1);
static int ghVec4INotEqual(GHvec v0, GHvec v1);

/* Vector Rounding and Clamping Operations */
static GHvec ghVecMin(GHvec v0, GHvec v1);
static GHvec ghVecMax(GHvec v0, GHvec v1);
static GHvec ghVecRound(GHvec v);
static GHvec ghVecTrunc(GHvec v);
static GHvec ghVecFloor(GHvec v);
static GHvec ghVecCeil(GHvec v);
static GHvec ghVecClamp(GHvec v, GHvec min, GHvec max);
static GHvec ghVecSaturate(GHvec v);

/* Integer Vector Bitwise Logical Operations */
static GHvec ghVecIAnd(GHvec v0, GHvec v1);
static GHvec ghVecIAndC(GHvec v0, GHvec v1);
static GHvec ghVecIOr(GHvec v0, GHvec v1);
static GHvec ghVecINor(GHvec v0, GHvec v1);
static GHvec ghVecIXor(GHvec v0, GHvec v1);

/* Vector Misc Operations */
static GHvec ghVecSelectMask(uint32_t idx0, uint32_t idx1, uint32_t idx2, uint32_t idx3);
static GHvec ghVecSelect(GHvec v0, GHvec v1, GHvec ctrl);

/* Scalar Computation Operations */
static float ghScalarSin(float rad);
static float ghScalarCos(float rad);
static void ghScalarSinCos(float *s, float *c, float rad);

/* Vector Computation Operations */
static GHvec ghVecAdd(GHvec v0, GHvec v1);
static GHvec ghVecSub(GHvec v0, GHvec v1);
static GHvec ghVecMul(GHvec v0, GHvec v1);
static GHvec ghVecMad(GHvec v0, GHvec v1, GHvec v2);
static GHvec ghVecDiv(GHvec v0, GHvec v1);
static GHvec ghVecScale(GHvec v, float s);
static GHvec ghVecRcp(GHvec v);
static GHvec ghVecRcpF(GHvec v);
static GHvec ghVecSqrt(GHvec v);
static GHvec ghVecRcpSqrt(GHvec v);
static GHvec ghVecRcpSqrtF(GHvec v);
static GHvec ghVecNeg(GHvec v);
static GHvec ghVecAbs(GHvec v);
static GHvec ghVecModRadians(GHvec v);
static GHvec ghVecSin(GHvec v);
static GHvec ghVecCos(GHvec v);
static void ghVecSinCos(GHvec *dstSin, GHvec *dstCos, GHvec v);

/* VectorN Computation Operations */
static GHvec ghVec2Dot(GHvec v0, GHvec v1);
static GHvec ghVec3Dot(GHvec v0, GHvec v1);
static GHvec ghVec4Dot(GHvec v0, GHvec v1);
static GHvec ghVec3Cross(GHvec v0, GHvec v1);
static GHvec ghVec2Length(GHvec v);
static GHvec ghVec3Length(GHvec v);
static GHvec ghVec4Length(GHvec v);
static GHvec ghVec2NormalizeF(GHvec v);
static GHvec ghVec3NormalizeF(GHvec v);
static GHvec ghVec4NormalizeF(GHvec v);
static GHvec ghVec4Transform(GHvec v, const GHmat *m);

/* Matrix Computation Operations */
static GHmat ghMatSet(float, float, float, float, float, float, float, float, float, float, float, float, float, float, float, float);
static GHmat ghMatIdentity(void);
static GHmat ghMatTranslation(float x, float y, float z);
static GHmat ghMatScaling(float x, float y, float z);
static GHmat ghMatRotationX(float rad);
static GHmat ghMatRotationY(float rad);
static GHmat ghMatRotationZ(float rad);
static GHmat ghMatRotationQuat(GHvec q);
static GHmat ghMatTranspose(const GHmat *m);
static GHmat ghMatMul(const GHmat *m0, const GHmat *m1);
static GHmat ghMatPerspective(float fovy, float aspect, float n, float f);
static GHmat ghMatLookAt(GHvec eyePos, GHvec lookPos, GHvec upDir);

/* Quaternion Computation Operations */
static GHvec ghQuatMul(GHvec q0, GHvec q1);

/*****************************************************************************
 *
 * Implementation
 *
 *****************************************************************************/

/* Helper Constants */
extern const GHvecF32 nGHOne;
extern const GHvecF32 nGHZero;
extern const GHvecS32 nGHAbsMask;
extern const GHvecF32 nGHNoFraction;
extern const GHvecF32 nGHPi;
extern const GHvecF32 nGHRcpPi;
extern const GHvecF32 nGHTwoPi;
extern const GHvecF32 nGHRcpTwoPi;
extern const GHvecF32 nGHSinCoefficients0;
extern const GHvecF32 nGHSinCoefficients1;
extern const GHvecF32 nGHSinCoefficients2;
extern const GHvecF32 nGHCosCoefficients0;
extern const GHvecF32 nGHCosCoefficients1;
extern const GHvecF32 nGHCosCoefficients2;
extern const GHvecF32 nGHIdentityR0;
extern const GHvecF32 nGHIdentityR1;
extern const GHvecF32 nGHIdentityR2;
extern const GHvecF32 nGHIdentityR3;
extern const GHvecF32 nGHNegateX;
extern const GHvecF32 nGHNegateY;
extern const GHvecF32 nGHNegateZ;
extern const GHvecF32 nGHNegateW;
extern const GHvecS32 nGHMask2;
extern const GHvecS32 nGHMask3;
extern const GHvecF32 nGHOne3;
extern const GHvecS32 nGHMaskX;
extern const GHvecS32 nGHMaskY;
extern const GHvecS32 nGHMaskZ;
extern const GHvecS32 nGHMaskW;
extern const GHvecF32 nGHNegOne;
extern const GHvecS32 nGHNegOneMask;
extern const GHvecS32 nGHOneHalfMinusEpsilon;
extern const GHvecF32 nGHEpsilon;

/* FloatN Set Operations */
static GH_INLINE GHfloat2 ghFloat2Set(float x, float y) {
    GHfloat2 res;
    res.x = x; res.y = y;
    return res;
}
static GH_INLINE GHfloat2 ghFloat2Setv(const float *v) {
    GHfloat2 res;
    res.x = v[0]; res.y = v[1];
    return res;
}
static GH_INLINE GHfloat2A ghFloat2ASet(float x, float y) {
    GHfloat2A res;
    res.x = x; res.y = y;
    return res;
}
static GH_INLINE GHfloat2A ghFloat2ASetv(const float *v) {
    GHfloat2A res;
    res.x = v[0]; res.y = v[1];
    return res;
}
static GH_INLINE GHfloat3 ghFloat3Set(float x, float y, float z) {
    GHfloat3 res;
    res.x = x; res.y = y; res.z = z;
    return res;
}
static GH_INLINE GHfloat3 ghFloat3Setv(const float *v) {
    GHfloat3 res;
    res.x = v[0]; res.y = v[1]; res.z = v[2];
    return res;
}
static GH_INLINE GHfloat3A ghFloat3ASet(float x, float y, float z) {
    GHfloat3A res;
    res.x = x; res.y = y; res.z = z;
    return res;
}
static GH_INLINE GHfloat3A ghFloat3ASetv(const float *v) {
    GHfloat3A res;
    res.x = v[0]; res.y = v[1]; res.z = v[2];
    return res;
}
static GH_INLINE GHfloat4 ghFloat4Set(float x, float y, float z, float w) {
    GHfloat4 res;
    res.x = x; res.y = y; res.z = z; res.w = w;
    return res;
}
static GH_INLINE GHfloat4 ghFloat4Setv(const float *v) {
    GHfloat4 res;
    res.x = v[0]; res.y = v[1]; res.z = v[2]; res.w = v[3];
    return res;
}
static GH_INLINE GHfloat4A ghFloat4ASet(float x, float y, float z, float w) {
    GHfloat4A res;
    res.x = x; res.y = y; res.z = z; res.w = w;
    return res;
}
static GH_INLINE GHfloat4A ghFloat4ASetv(const float *v) {
    GHfloat4A res;
    res.x = v[0]; res.y = v[1]; res.z = v[2]; res.w = v[3];
    return res;
}
/* FloatMxN Set Operations */
static GH_INLINE GHfloat3x3 ghFloat3x3Set(
        float m00, float m01, float m02,
        float m10, float m11, float m12,
        float m20, float m21, float m22) {
    GHfloat3x3 res;
    res._00 = m00; res._01 = m01; res._02 = m02;
    res._10 = m10; res._11 = m11; res._12 = m12;
    res._20 = m20; res._21 = m21; res._22 = m22;
    return res;
}
static GH_INLINE GHfloat3x3 ghFloat3x3Setv(const float *v) {
    GHfloat3x3 res;
    res._00 = v[0]; res._01 = v[1];  res._02 = v[2];
    res._10 = v[3]; res._11 = v[4];  res._12 = v[5];
    res._20 = v[6]; res._21 = v[7];  res._22 = v[8];
    return res;
}
static GH_INLINE GHfloat4x3 ghFloat4x3Set(
        float m00, float m01, float m02,
        float m10, float m11, float m12,
        float m20, float m21, float m22,
        float m30, float m31, float m32) {
    GHfloat4x3 res;
    res._00 = m00; res._01 = m01; res._02 = m02;
    res._10 = m10; res._11 = m11; res._12 = m12;
    res._20 = m20; res._21 = m21; res._22 = m22;
    res._30 = m30; res._31 = m31; res._32 = m32;
    return res;
}
static GH_INLINE GHfloat4x3 ghFloat4x3Setv(const float *v) {
    GHfloat4x3 res;
    res._00 = v[0]; res._01 = v[1];  res._02 = v[2];
    res._10 = v[3]; res._11 = v[4];  res._12 = v[5];
    res._20 = v[6]; res._21 = v[7];  res._22 = v[8];
    res._30 = v[9]; res._31 = v[10]; res._32 = v[11];
    return res;
}
static GH_INLINE GHfloat4x3A ghFloat4x3ASet(
        float m00, float m01, float m02,
        float m10, float m11, float m12,
        float m20, float m21, float m22,
        float m30, float m31, float m32) {
    GHfloat4x3A res;
    res._00 = m00; res._01 = m01; res._02 = m02;
    res._10 = m10; res._11 = m11; res._12 = m12;
    res._20 = m20; res._21 = m21; res._22 = m22;
    res._30 = m30; res._31 = m31; res._32 = m32;
    return res;
}
static GH_INLINE GHfloat4x3A ghFloat4x3ASetv(const float *v) {
    GHfloat4x3A res;
    res._00 = v[0];  res._01 = v[1];  res._02 = v[2];
    res._10 = v[3];  res._11 = v[4];  res._12 = v[5];
    res._20 = v[6];  res._21 = v[7];  res._22 = v[8];
    res._30 = v[9]; res._31 = v[10]; res._32 = v[11];
    return res;
}
static GH_INLINE GHfloat4x4 ghFloat4x4Set(
        float m00, float m01, float m02, float m03,
        float m10, float m11, float m12, float m13,
        float m20, float m21, float m22, float m23,
        float m30, float m31, float m32, float m33) {
    GHfloat4x4 res;
    res._00 = m00; res._01 = m01; res._02 = m02; res._03 = m03;
    res._10 = m10; res._11 = m11; res._12 = m12; res._13 = m13;
    res._20 = m20; res._21 = m21; res._22 = m22; res._23 = m23;
    res._30 = m30; res._31 = m31; res._32 = m32; res._33 = m33;
    return res;
}
static GH_INLINE GHfloat4x4 ghFloat4x4Setv(const float *v) {
    GHfloat4x4 res;
    res._00 = v[0];  res._01 = v[1];  res._02 = v[2];  res._03 = v[3];
    res._10 = v[4];  res._11 = v[5];  res._12 = v[6];  res._13 = v[7];
    res._20 = v[8];  res._21 = v[9];  res._22 = v[10]; res._23 = v[11];
    res._30 = v[12]; res._31 = v[13]; res._32 = v[14]; res._33 = v[15];
    return res;
}
static GH_INLINE GHfloat4x4A ghFloat4x4ASet(
        float m00, float m01, float m02, float m03,
        float m10, float m11, float m12, float m13,
        float m20, float m21, float m22, float m23,
        float m30, float m31, float m32, float m33) {
    GHfloat4x4A res;
    res._00 = m00; res._01 = m01; res._02 = m02; res._03 = m03;
    res._10 = m10; res._11 = m11; res._12 = m12; res._13 = m13;
    res._20 = m20; res._21 = m21; res._22 = m22; res._23 = m23;
    res._30 = m30; res._31 = m31; res._32 = m32; res._33 = m33;
    return res;
}
static GH_INLINE GHfloat4x4A ghFloat4x4ASetv(const float *v) {
    GHfloat4x4A res;
    res._00 = v[0];  res._01 = v[1];  res._02 = v[2];  res._03 = v[3];
    res._10 = v[4];  res._11 = v[5];  res._12 = v[6];  res._13 = v[7];
    res._20 = v[8];  res._21 = v[9];  res._22 = v[10]; res._23 = v[11];
    res._30 = v[12]; res._31 = v[13]; res._32 = v[14]; res._33 = v[15];
    return res;
}
/* Vector Load Operations */
static GH_INLINE GHvec ghVecLoadFloat(const float *src) {
    return _mm_load_ss(src);
}
static GH_INLINE GHvec ghVecLoadFloat2(const GHfloat2 *src) {
    GHvec x, y;
    x = _mm_load_ss(&src->x);
    y = _mm_load_ss(&src->y);
    return _mm_unpacklo_ps(x, y);
}
static GH_INLINE GHvec ghVecLoadFloat3(const GHfloat3 *src) {
    GHvec x = _mm_load_ss(&src->x);
    GHvec y = _mm_load_ss(&src->y);
    GHvec z = _mm_load_ss(&src->z);
    GHvec xy = _mm_unpacklo_ps(x, y);
    return _mm_movelh_ps(xy, z);
}
static GH_INLINE GHvec ghVecLoadFloat4(const GHfloat4 *src) {
    return _mm_loadu_ps(&src->x);
}
static GH_INLINE GHvec ghVecLoadFloat2A(const GHfloat2A *src) {
    __m128i v = _mm_loadl_epi64((const __m128i *) src);
    return *((const GHvec *) &v);
}
static GH_INLINE GHvec ghVecLoadFloat3A(const GHfloat3A *src) {
    return _mm_load_ps(&src->x);
}
static GH_INLINE GHvec ghVecLoadFloat4A(const GHfloat4A *src) {
    return _mm_load_ps(&src->x);
}
/* Matrix Load Operations */
static GH_INLINE GHmat ghMatLoadFloat3x3(const GHfloat3x3 *src) {
    GHmat m;
    GHvec v1, v2, v3, z, t1, t2, t3, t4, t5;

    z = _mm_setzero_ps();

    v1 = _mm_loadu_ps(&src->m[0][0]);
    v2 = _mm_loadu_ps(&src->m[1][1]);
    v3 = _mm_load_ss(&src->m[2][2]);

    t1 = _mm_unpackhi_ps(v1, z);
    t2 = _mm_unpacklo_ps(v2, z);
    t3 = _mm_shuffle_ps(v3, t2, _MM_SHUFFLE(0, 1, 0, 0));
    t4 = _mm_movehl_ps(t2, t3);
    t5 = _mm_movehl_ps(z, t1);

    m.r[0] = _mm_movelh_ps(v1, t1);
    m.r[1] = _mm_add_ps(t4, t5);
    m.r[2] = _mm_shuffle_ps(v2, v3, _MM_SHUFFLE(1, 0, 3, 2));
    m.r[3] = nGHIdentityR3.v;

    return m;
}
static GH_INLINE GHmat ghMatLoadFloat4x3(const GHfloat4x3 *src) {
    GHmat m;
    GHvec t0, t1, t2, t3;
    __m128i t3i;
    t0 = _mm_loadu_ps(&src->m[0][0]);
    t1 = _mm_loadu_ps(&src->m[1][1]);
    t3 = _mm_loadu_ps(&src->m[2][2]);
    t2 = _mm_shuffle_ps(t1, t3, _MM_SHUFFLE(0, 0, 3, 2));
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(3, 3, 1, 0));
    t1 = _mm_shuffle_ps(t1, t1, _MM_SHUFFLE(1, 1, 0, 2));
    t0 = _mm_and_ps(t0, nGHMask3.v);
    t1 = _mm_and_ps(t1, nGHMask3.v);
    t2 = _mm_and_ps(t2, nGHMask3.v);
    t3i = _mm_srli_si128(*((const __m128i *) &t3), 4);
    t3i = _mm_or_si128(t3i, *((const __m128i *) &nGHIdentityR3.v));
    m.r[0] = t0;
    m.r[1] = t1;
    m.r[2] = t2;
    m.r[3] = *((const GHvec *) &t3i);
    return m;
}
static GH_INLINE GHmat ghMatLoadFloat4x4(const GHfloat4x4 *src) {
    GHmat m;
    m.r[0] = _mm_loadu_ps(&src->_00);
    m.r[1] = _mm_loadu_ps(&src->_10);
    m.r[2] = _mm_loadu_ps(&src->_20);
    m.r[3] = _mm_loadu_ps(&src->_30);
    return m;
}
static GH_INLINE GHmat ghMatLoadFloat4x3A(const GHfloat4x3A *src) {
    GHmat m;
    GHvec t0, t1, t2, t3;
    __m128i t3i;
    t0 = _mm_load_ps(&src->m[0][0]);
    t1 = _mm_load_ps(&src->m[1][1]);
    t3 = _mm_load_ps(&src->m[2][2]);
    t2 = _mm_shuffle_ps(t1, t3, _MM_SHUFFLE(0, 0, 3, 2));
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(3, 3, 1, 0));
    t1 = _mm_shuffle_ps(t1, t1, _MM_SHUFFLE(1, 1, 0, 2));
    t0 = _mm_and_ps(t0, nGHMask3.v);
    t1 = _mm_and_ps(t1, nGHMask3.v);
    t2 = _mm_and_ps(t2, nGHMask3.v);
    t3i = _mm_srli_si128(*((const __m128i *) &t3), 4);
    t3i = _mm_or_si128(t3i, *((const __m128i *) &nGHIdentityR3.v));
    m.r[0] = t0;
    m.r[1] = t1;
    m.r[2] = t2;
    m.r[3] = *((const GHvec *) &t3i);
    return m;
}
static GH_INLINE GHmat ghMatLoadFloat4x4A(const GHfloat4x4A *src) {
    GHmat m;
    m.r[0] = _mm_load_ps(&src->_00);
    m.r[1] = _mm_load_ps(&src->_10);
    m.r[2] = _mm_load_ps(&src->_20);
    m.r[3] = _mm_load_ps(&src->_30);
    return m;
}
/* Vector Store Operations */
static GH_INLINE void ghVecStoreFloat(float *dst, GHvec v) {
    _mm_store_ss(dst, v);
}
static GH_INLINE void ghVecStoreFloat2(GHfloat2 *dst, GHvec v) {
    GHvec yyyy = _mm_shuffle_ps(v, v, _MM_SHUFFLE(1, 1, 1, 1));
    _mm_store_ss(&dst->x, v);
    _mm_store_ss(&dst->y, yyyy);
}
static GH_INLINE void ghVecStoreFloat3(GHfloat3 *dst, GHvec v) {
    GHvec yyyy = _mm_shuffle_ps(v, v, _MM_SHUFFLE(1, 1, 1, 1));
    GHvec zzzz = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 2, 2, 2));
    _mm_store_ss(&dst->x, v);
    _mm_store_ss(&dst->y, yyyy);
    _mm_store_ss(&dst->z, zzzz);
}
static GH_INLINE void ghVecStoreFloat4(GHfloat4 *dst, GHvec v) {
    _mm_storeu_ps(&dst->x, v);
}
static GH_INLINE void ghVecStoreFloat2A(GHfloat2A *dst, GHvec v) {
    _mm_storel_epi64((__m128i *) dst, *((const __m128i *) &v));
}
static GH_INLINE void ghVecStoreFloat3A(GHfloat3A *dst, GHvec v) {
    GHvec zzzz = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 2, 2, 2));
    _mm_storel_epi64((__m128i *) dst, *((const __m128i *) &v));
    _mm_store_ss(&dst->z, zzzz);
}
static GH_INLINE void ghVecStoreFloat4A(GHfloat4A *dst, GHvec v) {
    _mm_store_ps(&dst->x, v);
}
/* Matrix Store Operations */
static GH_INLINE void ghMatStoreFloat3x3(GHfloat3x3 *dst, const GHmat *m) {
    GHvec t0, t1, t2, work;
    t0 = m->r[0];
    t1 = m->r[1];
    t2 = m->r[2];
    work = _mm_shuffle_ps(t0, t1, _MM_SHUFFLE(0, 0, 2, 2));
    t0 = _mm_shuffle_ps(t0, work, _MM_SHUFFLE(2, 0, 1, 0));
    _mm_storeu_ps(&dst->m[0][0], t0);
    t1 = _mm_shuffle_ps(t1, t2, _MM_SHUFFLE(1, 0, 2, 1));
    _mm_storeu_ps(&dst->m[1][1],t1);
    t2 = _mm_shuffle_ps(t2, t2, _MM_SHUFFLE(2, 2, 2, 2));
    _mm_store_ss(&dst->m[2][2], t2);
}
static GH_INLINE void ghMatStoreFloat4x3(GHfloat4x3 *dst, const GHmat *m) {
    GHvec t0, t1, t2, t3, t2x;
    t0 = m->r[0];
    t1 = m->r[1];
    t2 = m->r[2];
    t3 = m->r[3];
    t2x = _mm_shuffle_ps(t1, t2, _MM_SHUFFLE(1, 0, 2, 1));
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(2, 2, 0, 0));
    t0 = _mm_shuffle_ps(t0, t1, _MM_SHUFFLE(0, 2, 1, 0));
    t2 = _mm_shuffle_ps(t2, t3, _MM_SHUFFLE(0, 0, 2, 2));
    t2 = _mm_shuffle_ps(t2, t3, _MM_SHUFFLE(2, 1, 2, 0));
    _mm_storeu_ps(&dst->m[0][0], t0);
    _mm_storeu_ps(&dst->m[1][1], t2x);
    _mm_storeu_ps(&dst->m[2][2], t2);
}
static GH_INLINE void ghMatStoreFloat4x4(GHfloat4x4 *dst, const GHmat *m) {
    _mm_storeu_ps(&dst->_00, m->r[0]);
	_mm_storeu_ps(&dst->_10, m->r[1]);
	_mm_storeu_ps(&dst->_20, m->r[2]);
	_mm_storeu_ps(&dst->_30, m->r[3]);
}
static GH_INLINE void ghMatStoreFloat4x3A(GHfloat4x3A *dst, const GHmat *m) {
    GHvec t, t0, t1, t2, t3;
    t0 = m->r[0];
    t1 = m->r[1];
    t2 = m->r[2];
    t3 = m->r[3];
    t = _mm_shuffle_ps(t0, t1, _MM_SHUFFLE(1, 0, 2, 2));
    t1 = _mm_shuffle_ps(t1, t2, _MM_SHUFFLE(1, 0, 2, 1));
    t0 = _mm_shuffle_ps(t0, t, _MM_SHUFFLE(2, 0, 1, 0));
    t2 = _mm_shuffle_ps(t2, t3, _MM_SHUFFLE(0, 0, 2, 2));
    t2 = _mm_shuffle_ps(t2, t3, _MM_SHUFFLE(2, 1, 2, 0));
    _mm_store_ps(&dst->m[0][0], t0);
    _mm_store_ps(&dst->m[1][1], t1);
    _mm_store_ps(&dst->m[2][2], t2);
}
static GH_INLINE void ghMatStoreFloat4x4A(GHfloat4x4A *dst, const GHmat *m) {
    _mm_store_ps(&dst->_00, m->r[0]);
    _mm_store_ps(&dst->_10, m->r[1]);
    _mm_store_ps(&dst->_20, m->r[2]);
    _mm_store_ps(&dst->_30, m->r[3]);
}
/* Vector Set and Get Operations */
static GH_INLINE GHvec ghVecZero(void) {
    return _mm_setzero_ps();
}
static GH_INLINE GHvec ghVecOne(void) {
    return nGHOne.v;
}
static GH_INLINE GHvec ghVecEpsilon(void) {
    return nGHEpsilon.v;
}
static GH_INLINE GHvec ghVecSet(float x, float y, float z, float w) {
    return _mm_set_ps(w, z, y, x);
}
static GH_INLINE GHvec ghVecReplicate(float s) {
    return _mm_set_ps1(s);
}
static GH_INLINE GHvec ghVecReplicatev(const float *src) {
    return _mm_load_ps1(src);
}
static GH_INLINE float ghVecGetX(GHvec v) {
    return _mm_cvtss_f32(v);
}
static GH_INLINE float ghVecGetY(GHvec v) {
    GHvec r = _mm_shuffle_ps(v, v, _MM_SHUFFLE(1, 1, 1, 1));
    return _mm_cvtss_f32(r);
}
static GH_INLINE float ghVecGetZ(GHvec v) {
    GHvec r = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 2, 2, 2));
    return _mm_cvtss_f32(r);
}
static GH_INLINE float ghVecGetW(GHvec v) {
    GHvec r = _mm_shuffle_ps(v, v, _MM_SHUFFLE(3, 3, 3, 3));
    return _mm_cvtss_f32(r);
}
/* Integer Vector Set Operations */
static GH_INLINE GHvec ghVecITrue(void) {
    __m128i v = _mm_set1_epi32(-1);
    return *((const GHvec *) &v);
}
static GH_INLINE GHvec ghVecIFalse(void) {
    return _mm_setzero_ps();
}
static GH_INLINE GHvec ghVecISet(uint32_t x, uint32_t y, uint32_t z, uint32_t w) {
    __m128i v = _mm_set_epi32(w, z, y, x);
    return *((const GHvec *) &v);
}
static GH_INLINE GHvec ghVecIReplicate(uint32_t s) {
    __m128i v = _mm_set1_epi32(s);
    return *((const GHvec *) &v);
}
static GH_INLINE GHvec ghVecIReplicatev(const uint32_t *src) {
    return _mm_load_ps1((const float *) src);
}
/* Vector Comparison Operations */
static GH_INLINE GHvec ghVecEqual(GHvec v0, GHvec v1) {
    return _mm_cmpeq_ps(v0, v1);
}
static GH_INLINE GHvec ghVecNearEqual(GHvec v0, GHvec v1, GHvec epsilon) {
    GHvec d, t;
    d = _mm_sub_ps(v0, v1);
    t = _mm_setzero_ps();
    t = _mm_sub_ps(t, d);
    t = _mm_max_ps(t, d);
    t = _mm_cmple_ps(t, epsilon);
    return t;
}
static GH_INLINE GHvec ghVecNotEqual(GHvec v0, GHvec v1) {
    return _mm_cmpneq_ps(v0, v1);
}
static GH_INLINE GHvec ghVecLess(GHvec v0, GHvec v1) {
    return _mm_cmplt_ps(v0, v1);
}
static GH_INLINE GHvec ghVecLessEqual(GHvec v0, GHvec v1) {
    return _mm_cmple_ps(v0, v1);
}
static GH_INLINE GHvec ghVecGreater(GHvec v0, GHvec v1) {
    return _mm_cmpgt_ps(v0, v1);
}
static GH_INLINE GHvec ghVecGreaterEqual(GHvec v0, GHvec v1) {
    return _mm_cmpge_ps(v0, v1);
}
/* VectorN Comparison Operations */
static GH_INLINE int ghVec2Equal(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpeq_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x3) == 0x3;
}
static GH_INLINE int ghVec3Equal(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpeq_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x7) == 0x7;
}
static GH_INLINE int ghVec4Equal(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpeq_ps(v0, v1);
    return _mm_movemask_ps(v) == 0xf;
}
static GH_INLINE int ghVec2NearEqual(GHvec v0, GHvec v1, GHvec epsilon) {
    GHvec d, t;
    d = _mm_sub_ps(v0, v1);
    t = _mm_setzero_ps();
    t = _mm_sub_ps(t, d);
    t = _mm_max_ps(t, d);
    t = _mm_cmple_ps(t, epsilon);
    return (_mm_movemask_ps(t) & 0x3) == 0x3;
}
static GH_INLINE int ghVec3NearEqual(GHvec v0, GHvec v1, GHvec epsilon) {
    GHvec d, t;
    d = _mm_sub_ps(v0, v1);
    t = _mm_setzero_ps();
    t = _mm_sub_ps(t, d);
    t = _mm_max_ps(t, d);
    t = _mm_cmple_ps(t, epsilon);
    return (_mm_movemask_ps(t) & 0x7) == 0x7;
}
static GH_INLINE int ghVec4NearEqual(GHvec v0, GHvec v1, GHvec epsilon) {
    GHvec d, t;
    d = _mm_sub_ps(v0, v1);
    t = _mm_setzero_ps();
    t = _mm_sub_ps(t, d);
    t = _mm_max_ps(t, d);
    t = _mm_cmple_ps(t, epsilon);
    return _mm_movemask_ps(t) == 0xf;
}
static GH_INLINE int ghVec2NotEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpeq_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x3) != 0x3;
}
static GH_INLINE int ghVec3NotEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpeq_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x7) != 0x7;
}
static GH_INLINE int ghVec4NotEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpeq_ps(v0, v1);
    return _mm_movemask_ps(v) != 0xf;
}
static GH_INLINE int ghVec2Less(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmplt_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x3) == 0x3;
}
static GH_INLINE int ghVec3Less(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmplt_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x7) == 0x7;
}
static GH_INLINE int ghVec4Less(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmplt_ps(v0, v1);
    return _mm_movemask_ps(v) == 0xf;
}
static GH_INLINE int ghVec2LessEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmple_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x3) == 0x3;
}
static GH_INLINE int ghVec3LessEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmple_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x7) == 0x7;
}
static GH_INLINE int ghVec4LessEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmple_ps(v0, v1);
    return _mm_movemask_ps(v) == 0xf;
}
static GH_INLINE int ghVec2Greater(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpgt_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x3) == 0x3;
}
static GH_INLINE int ghVec3Greater(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpgt_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x7) == 0x7;
}
static GH_INLINE int ghVec4Greater(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpgt_ps(v0, v1);
    return _mm_movemask_ps(v) == 0xf;
}
static GH_INLINE int ghVec2GreaterEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpge_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x3) == 0x3;
}
static GH_INLINE int ghVec3GreaterEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpge_ps(v0, v1);
    return (_mm_movemask_ps(v) & 0x7) == 0x7;
}
static GH_INLINE int ghVec4GreaterEqual(GHvec v0, GHvec v1) {
    GHvec v = _mm_cmpge_ps(v0, v1);
    return _mm_movemask_ps(v) == 0xf;
}
/* Integer Vector Comparison Operations */
static GH_INLINE GHvec ghVecIEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return *((const GHvec *) &v);
}
static GH_INLINE GHvec ghVecINotEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return _mm_xor_ps(*((const GHvec *) &v), nGHNegOneMask.v);
}
/* Integer VectorN Comparison Operations */
static GH_INLINE int ghVec2IEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return (_mm_movemask_ps(*((const GHvec *) &v)) & 0x3) == 0x3;
}
static GH_INLINE int ghVec3IEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return (_mm_movemask_ps(*((const GHvec *) &v)) & 0x7) == 0x7;
}
static GH_INLINE int ghVec4IEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return _mm_movemask_ps(*((const GHvec *) &v)) == 0xf;
}
static GH_INLINE int ghVec2INotEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return (_mm_movemask_ps(*((const GHvec *) &v)) & 0x3) != 0x3;
}
static GH_INLINE int ghVec3INotEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return (_mm_movemask_ps(*((const GHvec *) &v)) & 0x7) != 0x7;
}
static GH_INLINE int ghVec4INotEqual(GHvec v0, GHvec v1) {
    __m128i v = _mm_cmpeq_epi32(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return _mm_movemask_ps(*((const GHvec *) &v)) != 0xf;
}
/* Vector Rounding and Clamping Operations */
static GH_INLINE GHvec ghVecMin(GHvec v0, GHvec v1) {
    return _mm_min_ps(v0, v1);
}
static GH_INLINE GHvec ghVecMax(GHvec v0, GHvec v1) {
    return _mm_max_ps(v0, v1);
}
static GH_INLINE GHvec ghVecRound(GHvec v) {
    __m128i t, i;
    GHvec r;
    t = _mm_and_si128(*((const __m128i *) &v), *((const __m128i *) &nGHAbsMask.v));
    t = _mm_cmplt_epi32(t, *((const __m128i *) &nGHNoFraction.v));
    i = _mm_cvtps_epi32(v);
    r = _mm_cvtepi32_ps(i);
    r = _mm_and_ps(r, *((const GHvec *) &t));
    t = _mm_andnot_si128(t, *((const __m128i *) &v));
    r = _mm_or_ps(r, *((const GHvec *) &t));
    return r;
}
static GH_INLINE GHvec ghVecTrunc(GHvec v) {
    __m128i t, i;
    GHvec r;
    t = _mm_and_si128(*((const __m128i *) &v), *((const __m128i *) &nGHAbsMask.v));
    t = _mm_cmplt_epi32(t, *((const __m128i *) &nGHNoFraction.v));
    i = _mm_cvttps_epi32(v);
    r = _mm_cvtepi32_ps(i);
    r = _mm_and_ps(r, *((const GHvec *) &t));
    t = _mm_andnot_si128(t, *((const __m128i *) &v));
    r = _mm_or_ps(r, *((const GHvec *) &t));
    return r;
}
static GH_INLINE GHvec ghVecFloor(GHvec v) {
    GHvec r;
    __m128i i;
    r = _mm_sub_ps(v, nGHOneHalfMinusEpsilon.v);
    i = _mm_cvtps_epi32(r);
    r = _mm_cvtepi32_ps(i);
	return r;
}
static GH_INLINE GHvec ghVecCeil(GHvec v) {
    GHvec r;
    __m128i i;
    r = _mm_add_ps(v, nGHOneHalfMinusEpsilon.v);
    i = _mm_cvtps_epi32(r);
    r = _mm_cvtepi32_ps(i);
	return r;
}
static GH_INLINE GHvec ghVecClamp(GHvec v, GHvec min, GHvec max) {
    GHvec r;
    r = _mm_max_ps(v, min);
	r = _mm_min_ps(r, max);
	return r;
}
static GH_INLINE GHvec ghVecSaturate(GHvec v) {
    GHvec r;
    r = _mm_max_ps(v, nGHZero.v);
    r = _mm_min_ps(r, nGHOne.v);
    return r;
}
/* Integer Vector Bitwise Logical Operations */
static GH_INLINE GHvec ghVecIAnd(GHvec v0, GHvec v1) {
    __m128i r = _mm_and_si128(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return *((const GHvec *) &r);
}
static GH_INLINE GHvec ghVecIAndC(GHvec v0, GHvec v1) {
    __m128i r = _mm_andnot_si128(*((const __m128i *) &v1), *((const __m128i *) &v0));
    return *((const GHvec *) &r);
}
static GH_INLINE GHvec ghVecIOr(GHvec v0, GHvec v1) {
    __m128i r = _mm_or_si128(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return *((const GHvec *) &r);
}
static GH_INLINE GHvec ghVecINor(GHvec v0, GHvec v1) {
    __m128i r;
    r = _mm_or_si128(*((const __m128i *) &v0), *((const __m128i *) &v1));
    r = _mm_andnot_si128(r, *((const __m128i *) &nGHNegOneMask.v));
    return *((const GHvec *) &r);
}
static GH_INLINE GHvec ghVecIXor(GHvec v0, GHvec v1) {
    __m128i r = _mm_xor_si128(*((const __m128i *) &v0), *((const __m128i *) &v1));
    return *((const GHvec *) &r);
}
/* Vector Misc Operations */
static GH_INLINE GHvec ghVecSelectMask(uint32_t idx0, uint32_t idx1, uint32_t idx2, uint32_t idx3) {
    __m128i t;
    t = _mm_set_epi32(idx3, idx2, idx1, idx0);
    t = _mm_cmpgt_epi32(t, *((const __m128i *) &nGHZero.v));
	return *((const GHvec *) &t);
}
static GH_INLINE GHvec ghVecSelect(GHvec v0, GHvec v1, GHvec ctrl) {
    GHvec t0, t1;
    t0 = _mm_andnot_ps(ctrl, v0);
    t1 = _mm_and_ps(v1, ctrl);
    t1 = _mm_or_ps(t0, t1);
    return t1;
}
/* Scalar Computation Operations */
static GH_INLINE float ghScalarSin(float rad) {
    return sinf(rad);
}
static GH_INLINE float ghScalarCos(float rad) {
    return cosf(rad);
}
static GH_INLINE void ghScalarSinCos(float *s, float *c, float rad) {
    *s = sinf(rad);
    *c = cosf(rad);
}
/* Vector Computation Operations */
static GH_INLINE GHvec ghVecAdd(GHvec v0, GHvec v1) {
    return _mm_add_ps(v0, v1);
}
static GH_INLINE GHvec ghVecSub(GHvec v0, GHvec v1) {
    return _mm_sub_ps(v0, v1);
}
static GH_INLINE GHvec ghVecMul(GHvec v0, GHvec v1) {
    return _mm_mul_ps(v0, v1);
}
static GH_INLINE GHvec ghVecMad(GHvec v0, GHvec v1, GHvec v2) {
    GHvec r = _mm_mul_ps(v0, v1);
    return _mm_add_ps(r, v2);
}
static GH_INLINE GHvec ghVecDiv(GHvec v0, GHvec v1) {
    return _mm_div_ps(v0, v1);
}
static GH_INLINE GHvec ghVecScale(GHvec v, float s) {
    GHvec r = _mm_set_ps1(s);
    return _mm_mul_ps(v, r);
}
static GH_INLINE GHvec ghVecRcp(GHvec v) {
    return _mm_div_ps(nGHOne.v, v);
}
static GH_INLINE GHvec ghVecRcpF(GHvec v) {
    return _mm_rcp_ps(v);
}
static GH_INLINE GHvec ghVecSqrt(GHvec v) {
    return _mm_sqrt_ps(v);
}
static GH_INLINE GHvec ghVecRcpSqrt(GHvec v) {
    GHvec r = _mm_sqrt_ps(v);
    return _mm_div_ps(nGHOne.v, r);
}
static GH_INLINE GHvec ghVecRcpSqrtF(GHvec v) {
    return _mm_rsqrt_ps(v);
}
static GH_INLINE GHvec ghVecNeg(GHvec v) {
    GHvec z = _mm_setzero_ps();
    return _mm_sub_ps(z, v);
}
static GH_INLINE GHvec ghVecAbs(GHvec v) {
    GHvec r = _mm_setzero_ps();
    r = _mm_sub_ps(r, v);
    return _mm_max_ps(r, v);
}
static GH_INLINE GHvec ghVecModRadians(GHvec v) {
    GHvec r = _mm_mul_ps(v, nGHRcpTwoPi.v);
    r = ghVecRound(r);
    r = _mm_mul_ps(r, nGHTwoPi.v);
    return _mm_sub_ps(v, r);
}
static GH_INLINE GHvec ghVecSin(GHvec v) {
    GHvec v2, power, constants, r;

    r = ghVecModRadians(v);
    v2  = _mm_mul_ps(r, r);

    /* V1^3 */
    power = _mm_mul_ps(r, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients0.f[1]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^5 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients0.f[2]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^7 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients0.f[3]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^9 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients1.f[0]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^11 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients1.f[1]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^13 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients1.f[2]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^15 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients1.f[3]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^17 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients2.f[0]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^19 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients2.f[1]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^21 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients2.f[2]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^23 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHSinCoefficients2.f[3]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    return r;
}
static GH_INLINE GHvec ghVecCos(GHvec v) {
    GHvec v2, power, constants, r;

    r = ghVecModRadians(v);
    v2 = _mm_mul_ps(r, r);

    /* V^2 */
    constants = _mm_load_ps1(&nGHCosCoefficients0.f[1]);
    constants = _mm_mul_ps(constants, v2);
    r = _mm_add_ps(constants, nGHOne.v);

    /* V^4 */
    power = _mm_mul_ps(v2, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients0.f[2]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^6 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients0.f[3]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^8 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients1.f[0]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^10 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients1.f[1]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^12 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients1.f[2]);
    constants = _mm_mul_ps(constants,power);
    r = _mm_add_ps(r, constants);

    /* V^14 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients1.f[3]);
    constants = _mm_mul_ps(constants,power);
    r = _mm_add_ps(r, constants);

    /* V^16 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients2.f[0]);
    constants = _mm_mul_ps(constants,power);
    r = _mm_add_ps(r, constants);

    /* V^18 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients2.f[1]);
    constants = _mm_mul_ps(constants,power);
    r = _mm_add_ps(r, constants);

    /* V^20 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients2.f[2]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    /* V^22 */
    power = _mm_mul_ps(power, v2);
    constants = _mm_load_ps1(&nGHCosCoefficients2.f[3]);
    constants = _mm_mul_ps(constants, power);
    r = _mm_add_ps(r, constants);

    return r;
}
static GH_INLINE void ghVecSinCos(GHvec *dstSin, GHvec *dstCos, GHvec v) {
    GHvec v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13;
    GHvec v14, v15, v16, v17, v18, v19, v20, v21, v22, v23;
    GHvec s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11;
    GHvec c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11;
    GHvec sinv, cosv;

    v1 = ghVecModRadians(v);

    /*
    sin(V) ~= V - V^3 / 3! + V^5 / 5! - V^7 / 7! + V^9 / 9! - V^11 / 11! + V^13 / 13! -
              V^15 / 15! + V^17 / 17! - V^19 / 19! + V^21 / 21! - V^23 / 23! (for -PI <= V < PI)
    cos(V) ~= 1 - V^2 / 2! + V^4 / 4! - V^6 / 6! + V^8 / 8! - V^10 / 10! + V^12 / 12! -
              V^14 / 14! + V^16 / 16! - V^18 / 18! + V^20 / 20! - V^22 / 22! (for -PI <= V < PI)
    */

    v2 = _mm_mul_ps(v1, v1);
    v3 = _mm_mul_ps(v2, v1);
    v4 = _mm_mul_ps(v2, v2);
    v5 = _mm_mul_ps(v3, v2);
    v6 = _mm_mul_ps(v3, v3);
    v7 = _mm_mul_ps(v4, v3);
    v8 = _mm_mul_ps(v4, v4);
    v9 = _mm_mul_ps(v5, v4);
    v10 = _mm_mul_ps(v5, v5);
    v11 = _mm_mul_ps(v6, v5);
    v12 = _mm_mul_ps(v6, v6);
    v13 = _mm_mul_ps(v7, v6);
    v14 = _mm_mul_ps(v7, v7);
    v15 = _mm_mul_ps(v8, v7);
    v16 = _mm_mul_ps(v8, v8);
    v17 = _mm_mul_ps(v9, v8);
    v18 = _mm_mul_ps(v9, v9);
    v19 = _mm_mul_ps(v10, v9);
    v20 = _mm_mul_ps(v10, v10);
    v21 = _mm_mul_ps(v11, v10);
    v22 = _mm_mul_ps(v11, v11);
    v23 = _mm_mul_ps(v12, v11);

    s1 = _mm_load_ps1(&nGHSinCoefficients0.f[1]);
    s2 = _mm_load_ps1(&nGHSinCoefficients0.f[2]);
    s3 = _mm_load_ps1(&nGHSinCoefficients0.f[3]);
    s4 = _mm_load_ps1(&nGHSinCoefficients1.f[0]);
    s5 = _mm_load_ps1(&nGHSinCoefficients1.f[1]);
    s6 = _mm_load_ps1(&nGHSinCoefficients1.f[2]);
    s7 = _mm_load_ps1(&nGHSinCoefficients1.f[3]);
    s8 = _mm_load_ps1(&nGHSinCoefficients2.f[0]);
    s9 = _mm_load_ps1(&nGHSinCoefficients2.f[1]);
    s10 = _mm_load_ps1(&nGHSinCoefficients2.f[2]);
    s11 = _mm_load_ps1(&nGHSinCoefficients2.f[3]);

    c1 = _mm_load_ps1(&nGHCosCoefficients0.f[1]);
    c2 = _mm_load_ps1(&nGHCosCoefficients0.f[2]);
    c3 = _mm_load_ps1(&nGHCosCoefficients0.f[3]);
    c4 = _mm_load_ps1(&nGHCosCoefficients1.f[0]);
    c5 = _mm_load_ps1(&nGHCosCoefficients1.f[1]);
    c6 = _mm_load_ps1(&nGHCosCoefficients1.f[2]);
    c7 = _mm_load_ps1(&nGHCosCoefficients1.f[3]);
    c8 = _mm_load_ps1(&nGHCosCoefficients2.f[0]);
    c9 = _mm_load_ps1(&nGHCosCoefficients2.f[1]);
    c10 = _mm_load_ps1(&nGHCosCoefficients2.f[2]);
    c11 = _mm_load_ps1(&nGHCosCoefficients2.f[3]);

    s1 = _mm_mul_ps(s1, v3);
    sinv = _mm_add_ps(s1, v1);
    sinv = ghVecMad(s2, v5, sinv);
    sinv = ghVecMad(s3, v7, sinv);
    sinv = ghVecMad(s4, v9, sinv);
    sinv = ghVecMad(s5, v11, sinv);
    sinv = ghVecMad(s6, v13, sinv);
    sinv = ghVecMad(s7, v15, sinv);
    sinv = ghVecMad(s8, v17, sinv);
    sinv = ghVecMad(s9, v19, sinv);
    sinv = ghVecMad(s10, v21, sinv);
    sinv = ghVecMad(s11, v23, sinv);

    cosv = _mm_mul_ps(c1, v2);
    cosv = _mm_add_ps(cosv, nGHOne.v);
    cosv = ghVecMad(c2, v4, cosv);
    cosv = ghVecMad(c3, v6, cosv);
    cosv = ghVecMad(c4, v8, cosv);
    cosv = ghVecMad(c5, v10, cosv);
    cosv = ghVecMad(c6, v12, cosv);
    cosv = ghVecMad(c7, v14, cosv);
    cosv = ghVecMad(c8, v16, cosv);
    cosv = ghVecMad(c9, v18, cosv);
    cosv = ghVecMad(c10, v20, cosv);
    cosv = ghVecMad(c11, v22, cosv);

    *dstSin = sinv;
    *dstCos = cosv;
}
/* VectorN Computation Operations */
static GH_INLINE GHvec ghVec2Dot(GHvec v0, GHvec v1) {
    GHvec t0, t1;
    t0 = _mm_mul_ps(v0, v1);
    t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(1, 1, 1, 1));
    t0 = _mm_add_ss(t0, t1);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(0, 0, 0, 0));
    return t0;
}
static GH_INLINE GHvec ghVec3Dot(GHvec v0, GHvec v1) {
    GHvec t0, t1;
    t0 = _mm_mul_ps(v0, v1);
    t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 1, 2, 1));
    t0 = _mm_add_ss(t0, t1);
    t1 = _mm_shuffle_ps(t1, t1, _MM_SHUFFLE(1, 1, 1, 1));
    t0 = _mm_add_ss(t0, t1);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(0, 0, 0, 0));
    return t0;
}
static GH_INLINE GHvec ghVec4Dot(GHvec v0, GHvec v1) {
    GHvec t0, t1 = v1;
    /* t0 = | w0 * w1 | z0 * z1 | y0 * y1 | x0 * x1 | */
    t0 = _mm_mul_ps(v0, v1);
    /* t1 = | y0 * y1 | x0 * x1 | --- | --- | */
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(1, 0, 0, 0));
    /* t0 = | y0 * y1 + w0 * w1 | x0 * x1 + z0 * z1 | --- | --- | */
    t0 = _mm_add_ps(t0, t1);
    /* t1 = | --- | y0 * y1 + w0 * w1 | --- | --- | */
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(0, 3, 0, 0));
    /* t0 = | --- | x0 * x1 + y0 * y1 + z0 * z1 + w0 * w1 | --- | --- | */
    t0 = _mm_add_ps(t0, t1);
    /* t0 = | dot | dot | dot | dot | */
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 2, 2, 2));
    return t0;
}
static GH_INLINE GHvec ghVec3Cross(GHvec v0, GHvec v1) {
    GHvec t0, t1, r;
    /* t0 = | w0 | x0 | z0 | y0 | */
    t0 = _mm_shuffle_ps(v0, v0, _MM_SHUFFLE(3, 0, 2, 1));
    /* t1 = | w1 | y1 | x1 | z1 | */
    t1 = _mm_shuffle_ps(v1, v1, _MM_SHUFFLE(3, 1, 0, 2));
    /* r = | w0 * w1 | x0 * y1 | z0 * x1 | y0 * z1 | */
    r = _mm_mul_ps(t0, t1);
    /* t0 = | w0 | y0 | x0 | z0 | */
    t0 = _mm_shuffle_ps(v0, v0, _MM_SHUFFLE(3, 1, 0, 2));
    /* t1 = | w1 | x1 | z1 | y1 | */
    t1 = _mm_shuffle_ps(v1, v1, _MM_SHUFFLE(3, 0, 2, 1));
    /* t0 = | w0 * w1 | y0 * x1 | x0 * z1 | z0 * y1 | */
    t0 = _mm_mul_ps(t0, t1);
    /* r = | w0 * w1 - w0 * w1 | x0 * y1 - y0 * x1 | z0 * x1 - x0 * z1 | y0 * z1 - z0 * y1 | */
    return _mm_sub_ps(r, t0);
}
static GH_INLINE GHvec ghVec2Length(GHvec v) {
    GHvec t0, t1;
    t0 = _mm_mul_ps(v, v);
    t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(1, 1, 1, 1));
    t0 = _mm_add_ss(t0, t1);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(0, 0, 0, 0));
    t0 = _mm_sqrt_ps(t0);
    return t0;
}
static GH_INLINE GHvec ghVec3Length(GHvec v) {
    GHvec t0, t1;
    t0 = _mm_mul_ps(v, v);
    t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 1, 2, 1));
    t0 = _mm_add_ss(t0, t1);
    t1 = _mm_shuffle_ps(t1, t1, _MM_SHUFFLE(1, 1, 1, 1));
    t0 = _mm_add_ss(t0, t1);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(0, 0, 0, 0));
    t0 = _mm_sqrt_ps(t0);
    return t0;
}
static GH_INLINE GHvec ghVec4Length(GHvec v) {
    GHvec t0, t1 = v;
    t0 = _mm_mul_ps(v, v);
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(1, 0, 0, 0));
    t0 = _mm_add_ps(t0, t1);
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(0, 3, 0, 0));
    t0 = _mm_add_ps(t0, t1);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 2, 2, 2));
    t0 = _mm_sqrt_ps(t0);
    return t0;
}
static GH_INLINE GHvec ghVec2NormalizeF(GHvec v) {
    GHvec t0, t1;
    t0 = _mm_mul_ps(v, v);
    t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(1, 1, 1, 1));
    t0 = _mm_add_ss(t0, t1);
    t0 = _mm_rsqrt_ss(t0);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(0, 0, 0, 0));
    t0 = _mm_mul_ps(v, t0);
    return t0;
}
static GH_INLINE GHvec ghVec3NormalizeF(GHvec v) {
    GHvec t0, t1;
    t0 = _mm_mul_ps(v, v);
    t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 1, 2, 1));
    t0 = _mm_add_ss(t0, t1);
    t1 = _mm_shuffle_ps(t1, t1, _MM_SHUFFLE(1, 1, 1, 1));
    t0 = _mm_add_ss(t0, t1);
    t0 = _mm_rsqrt_ss(t0);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(0, 0, 0, 0));
    t0 = _mm_mul_ps(v, t0);
    return t0;
}
static GH_INLINE GHvec ghVec4NormalizeF(GHvec v) {
    GHvec t0, t1 = v;
    t0 = _mm_mul_ps(v, v);
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(1, 0, 0, 0));
    t0 = _mm_add_ps(t0, t1);
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(0, 3, 0, 0));
    t0 = _mm_add_ps(t0, t1);
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 2, 2, 2));
    t0 = _mm_rsqrt_ps(t0);
    t0 = _mm_mul_ps(v, t0);
    return t0;
}
static GH_INLINE GHvec ghVec4Transform(GHvec v, const GHmat *m) {
    GHvec tX, tY, tZ, tW;
    tX = _mm_shuffle_ps(v, v, _MM_SHUFFLE(0, 0, 0, 0));
    tY = _mm_shuffle_ps(v, v, _MM_SHUFFLE(1, 1, 1, 1));
    tZ = _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 2, 2, 2));
    tW = _mm_shuffle_ps(v, v, _MM_SHUFFLE(3, 3, 3, 3));
    tX = _mm_mul_ps(tX, m->r[0]);
    tY = _mm_mul_ps(tY, m->r[1]);
    tZ = _mm_mul_ps(tZ, m->r[2]);
    tW = _mm_mul_ps(tW, m->r[3]);
    tX = _mm_add_ps(tX, tY);
    tZ = _mm_add_ps(tZ, tW);
    tX = _mm_add_ps(tX, tZ);
    return tX;
}
/* Matrix Computation Operations */
static GH_INLINE GHmat ghMatSet(
        float m00, float m01, float m02, float m03,
        float m10, float m11, float m12, float m13,
        float m20, float m21, float m22, float m23,
        float m30, float m31, float m32, float m33) {
    GHmat m;
    m.r[0] = ghVecSet(m00, m01, m02, m03);
    m.r[1] = ghVecSet(m10, m11, m12, m13);
    m.r[2] = ghVecSet(m20, m21, m22, m23);
    m.r[3] = ghVecSet(m30, m31, m32, m33);
    return m;
}
static GH_INLINE GHmat ghMatIdentity(void) {
    GHmat m;
    m.r[0] = nGHIdentityR0.v;
    m.r[1] = nGHIdentityR1.v;
    m.r[2] = nGHIdentityR2.v;
    m.r[3] = nGHIdentityR3.v;
    return m;
}
static GH_INLINE GHmat ghMatTranslation(float x, float y, float z) {
    /* 1 0 0 0
       0 1 0 0
       0 0 1 0
       x y z 1 */
    GHmat m;
    m.r[0] = nGHIdentityR0.v;
    m.r[1] = nGHIdentityR1.v;
    m.r[2] = nGHIdentityR2.v;
    m.r[3] = _mm_set_ps(1.f, z, y, x);
    return m;
}
static GH_INLINE GHmat ghMatScaling(float x, float y, float z) {
    /* x 0 0 0
       0 y 0 0
       0 0 z 0
       0 0 0 1 */
    GHmat m;
    m.r[0] = _mm_set_ps(0.f, 0.f, 0.f, x);
    m.r[1] = _mm_set_ps(0.f, 0.f, y, .0f);
    m.r[2] = _mm_set_ps(0.f, z, 0.f, .0f);
    m.r[3] = nGHIdentityR3.v;
    return m;
}
static GH_INLINE GHmat ghMatRotationX(float rad) {
    /* 1  0 0 0
       0  c s 0
       0 -s c 0
       0  0 0 1 */
    GHvec t0, t1;
    GHmat m;

    m.r[0] = nGHIdentityR0.v;

    /* t0 = | 0 | 0 | 0 | sin(rad) | */
    t0 = _mm_set_ss(sinf(rad));
    /* t1 = | 0 | 0 | 0 | cos(rad) | */
    t1 = _mm_set_ss(cosf(rad));

    /* t0 = | 0 | sin(rad) | cos(rad) | 0 | */
    t0 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(3, 0, 0, 3));
    m.r[1] = t0;

    /* t0 = | 0 | cos(rad) | sin(rad) | 0 | */
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(3, 1, 2, 0));
    /* t0 = | 0 | cos(rad) | -sin(rad) | 0 | */
    t0 = _mm_mul_ps(t0, nGHNegateY.v);
    m.r[2] = t0;
    m.r[3] = nGHIdentityR3.v;

    return m;
}
static GH_INLINE GHmat ghMatRotationY(float rad) {
    /* c 0 -s 0
       0 1  0 0
       s 0  c 0
       0 0  0 1 */
    GHvec t0, t1;
    GHmat m;

    /* t0 = | 0 | 0 | 0 | sin(rad) | */
    t0 = _mm_set_ss(sinf(rad));
    /* t1 = | 0 | 0 | 0 | cos(rad) | */
    t1 = _mm_set_ss(cosf(rad));

    /* t1 = | 0 | sin(rad) | 0 | cos(rad) | */
    t1 = _mm_shuffle_ps(t1, t0, _MM_SHUFFLE(3, 0, 3, 0));
    /* t1 = | 0 | -sin(rad) | 0 | cos(rad) | */
    t1 = _mm_mul_ps(t1, nGHNegateZ.v);
    m.r[0] = t1;
    m.r[1] = nGHIdentityR1.v;

    /* t0 = | 0 | cos(rad) | 0 | sin(rad) | */
    t0 = _mm_shuffle_ps(t0, t1, _MM_SHUFFLE(3, 0, 3, 0));
    m.r[2] = t0;
    m.r[3] = nGHIdentityR3.v;

    return m;
}
static GH_INLINE GHmat ghMatRotationZ(float rad) {
    /*  c s 0 0
       -s c 0 0
        0 0 1 0
        0 0 0 1 */
    GHvec t0, t1;
    GHmat m;

    /* t0 = | 0 | 0 | 0 | sin(rad) | */
    t0 = _mm_set_ss(sinf(rad));
    /* t1 = | 0 | 0 | 0 | cos(rad) | */
    t1 = _mm_set_ss(cosf(rad));

    /* t0 = | 0 | 0 | sin(rad) | cos(rad) | */
    t0 = _mm_unpacklo_ps(t1, t0);
    m.r[0] = t0;

    /* t0 = | 0 | 0 | cos(rad) | sin(rad) | */
    t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(3, 2, 0, 1));
    /* t0 = | 0 | 0 | cos(rad) | -sin(rad) | */
    t0 = _mm_mul_ps(t0, nGHNegateX.v);
    m.r[1] = t0;
    m.r[2] = nGHIdentityR2.v;
    m.r[3] = nGHIdentityR3.v;

    return m;
}
static GH_INLINE GHmat ghMatRotationQuat(GHvec q) {
    /*
    |  1 - 2yy - 2zz  |      2xy - 2zw  |      2xz + 2yw  |  0  |
    |      2xy + 2zw  |  1 - 2xx - 2zz  |      2yz - 2xw  |  0  |
    |      2xz - 2yw  |      2yz + 2xw  |  1 - 2xx - 2yy  |  0  |
    |              0  |              0  |              0  |  1  |
    */
    GHmat m;
    GHvec q0, q1, v0, v1, v2, r0, r1, r2;
    q0 = _mm_add_ps(q, q);
    q1 = _mm_mul_ps(q, q0);
    v0 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(3, 0, 0, 1));
    v0 = _mm_and_ps(v0, nGHMask3.v);
    v1 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(3, 1, 2, 2));
    v1 = _mm_and_ps(v1, nGHMask3.v);
    r0 = _mm_sub_ps(nGHOne3.v, v0);
    r0 = _mm_sub_ps(r0, v1);
    v0 = _mm_shuffle_ps(q, q, _MM_SHUFFLE(3, 1, 0, 0));
    v1 = _mm_shuffle_ps(q0, q0, _MM_SHUFFLE(3, 2, 1, 2));
    v0 = _mm_mul_ps(v0, v1);
    v1 = _mm_shuffle_ps(q, q, _MM_SHUFFLE(3, 3, 3, 3));
    v2 = _mm_shuffle_ps(q0, q0, _MM_SHUFFLE(3, 0, 2, 1));
    v1 = _mm_mul_ps(v1, v2);
    r1 = _mm_sub_ps(v0, v1);
    r2 = _mm_add_ps(v0, v1);
    v0 = _mm_shuffle_ps(r1, r2, _MM_SHUFFLE(1, 2, 0, 1));
    v0 = _mm_shuffle_ps(v0, v0, _MM_SHUFFLE(3, 1, 0, 2));
    v1 = _mm_shuffle_ps(r1, r2, _MM_SHUFFLE(0, 0, 2, 2));
    v1 = _mm_shuffle_ps(v1, v1, _MM_SHUFFLE(3, 1, 2, 0));
    q1 = _mm_shuffle_ps(v0, r0, _MM_SHUFFLE(3, 0, 3, 2));
    q1 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(3, 0, 1, 2));
    m.r[0] = q1;
    q1 = _mm_shuffle_ps(v0, r0, _MM_SHUFFLE(3, 1, 1, 0));
    q1 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(3, 0, 2, 1));
    m.r[1] = q1;
    q1 = _mm_shuffle_ps(v1, r0, _MM_SHUFFLE(3, 2, 0, 1));
    m.r[2] = q1;
    m.r[3] = nGHIdentityR3.v;
    return m;
}
static GH_INLINE GHmat ghMatTranspose(const GHmat *m) {
    GHvec t0, t1, t2, t3;
    GHmat res;

    t0 = _mm_shuffle_ps(m->r[0], m->r[1], _MM_SHUFFLE(1, 0, 1, 0));
    t2 = _mm_shuffle_ps(m->r[0], m->r[1], _MM_SHUFFLE(3, 2, 3, 2));
    t1 = _mm_shuffle_ps(m->r[2], m->r[3], _MM_SHUFFLE(1, 0, 1, 0));
    t3 = _mm_shuffle_ps(m->r[2], m->r[3], _MM_SHUFFLE(3, 2, 3, 2));

    res.r[0] = _mm_shuffle_ps(t0, t1, _MM_SHUFFLE(2, 0, 2, 0));
    res.r[1] = _mm_shuffle_ps(t0, t1, _MM_SHUFFLE(3, 1, 3, 1));
    res.r[2] = _mm_shuffle_ps(t2, t3, _MM_SHUFFLE(2, 0, 2, 0));
    res.r[3] = _mm_shuffle_ps(t2, t3, _MM_SHUFFLE(3, 1, 3, 1));

	return res;
}
static GH_INLINE GHmat ghMatLookAt(GHvec eyePos, GHvec lookPos, GHvec upDir) {
    GHmat m;
    GHvec axisX, axisY, axisZ, negEyePos, d0, d1, d2;

    axisZ = _mm_sub_ps(eyePos, lookPos);
    axisZ = ghVec3NormalizeF(axisZ);
    axisX = ghVec3Cross(upDir, axisZ);
    axisY = ghVec3Cross(axisZ, axisX);

    negEyePos = _mm_mul_ps(nGHNegOne.v, eyePos);
    d0 = ghVec3Dot(negEyePos, axisX);
    d1 = ghVec3Dot(negEyePos, axisY);
    d2 = ghVec3Dot(negEyePos, axisZ);

    axisX = _mm_and_ps(axisX, nGHMask3.v);
    axisY = _mm_and_ps(axisY, nGHMask3.v);
    axisZ = _mm_and_ps(axisZ, nGHMask3.v);
    d0 = _mm_and_ps(d0, nGHMaskW.v);
    d1 = _mm_and_ps(d1, nGHMaskW.v);
    d2 = _mm_and_ps(d2, nGHMaskW.v);
    axisX = _mm_or_ps(axisX, d0);
    axisY = _mm_or_ps(axisY, d1);
    axisZ = _mm_or_ps(axisZ, d2);

    m.r[0] = axisX;
    m.r[1] = axisY;
    m.r[2] = axisZ;
    m.r[3] = nGHIdentityR3.v;
    return ghMatTranspose(&m);
}
static GH_INLINE GHmat ghMatPerspective(float fovy, float aspect, float n, float f) {
    GHmat m;
    GHvec values, t;
    float hfovy = fovy * .5f;
    float sy = cosf(hfovy) / sinf(hfovy);
    float range = f - n;
    values = ghVecSet(sy / aspect, sy, -(f+n)/range, -2.f*f*n/range);

    t = _mm_setzero_ps();
    t = _mm_move_ss(t, values);
    m.r[0] = t;

    t = values;
    t = _mm_and_ps(t, nGHMaskY.v);
    m.r[1] = t;

    t = _mm_shuffle_ps(values, values, _MM_SHUFFLE(2, 2, 2, 2));
    t = _mm_and_ps(t, nGHMaskZ.v);
    t = _mm_sub_ps(t, nGHIdentityR3.v);
    m.r[2] = t;

    t = _mm_shuffle_ps(values, values, _MM_SHUFFLE(3, 3, 3, 3));
    t = _mm_and_ps(t, nGHMaskZ.v);
    m.r[3] = t;

    return m;
}
static GH_INLINE GHmat ghMatMul(const GHmat *m0, const GHmat *m1) {
    GHmat res;
    GHvec tX, tY, tZ, tW;

    tW = m0->r[0];
    tX = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(0, 0, 0, 0));
    tY = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(1, 1, 1, 1));
    tZ = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(2, 2, 2, 2));
    tW = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(3, 3, 3, 3));
    tX = _mm_mul_ps(tX, m1->r[0]);
    tY = _mm_mul_ps(tY, m1->r[1]);
    tZ = _mm_mul_ps(tZ, m1->r[2]);
    tW = _mm_mul_ps(tW, m1->r[3]);
    tX = _mm_add_ps(tX, tZ);
    tY = _mm_add_ps(tY, tW);
    tX = _mm_add_ps(tX, tY);
    res.r[0] = tX;

    tW = m0->r[1];
    tX = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(0, 0, 0, 0));
    tY = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(1, 1, 1, 1));
    tZ = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(2, 2, 2, 2));
    tW = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(3, 3, 3, 3));
    tX = _mm_mul_ps(tX, m1->r[0]);
    tY = _mm_mul_ps(tY, m1->r[1]);
    tZ = _mm_mul_ps(tZ, m1->r[2]);
    tW = _mm_mul_ps(tW, m1->r[3]);
    tX = _mm_add_ps(tX, tZ);
    tY = _mm_add_ps(tY, tW);
    tX = _mm_add_ps(tX, tY);
    res.r[1] = tX;

    tW = m0->r[2];
    tX = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(0, 0, 0, 0));
    tY = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(1, 1, 1, 1));
    tZ = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(2, 2, 2, 2));
    tW = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(3, 3, 3, 3));
    tX = _mm_mul_ps(tX, m1->r[0]);
    tY = _mm_mul_ps(tY, m1->r[1]);
    tZ = _mm_mul_ps(tZ, m1->r[2]);
    tW = _mm_mul_ps(tW, m1->r[3]);
    tX = _mm_add_ps(tX, tZ);
    tY = _mm_add_ps(tY, tW);
    tX = _mm_add_ps(tX, tY);
    res.r[2] = tX;

    tW = m0->r[3];
    tX = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(0, 0, 0, 0));
    tY = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(1, 1, 1, 1));
    tZ = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(2, 2, 2, 2));
    tW = _mm_shuffle_ps(tW, tW, _MM_SHUFFLE(3, 3, 3, 3));
    tX = _mm_mul_ps(tX, m1->r[0]);
    tY = _mm_mul_ps(tY, m1->r[1]);
    tZ = _mm_mul_ps(tZ, m1->r[2]);
    tW = _mm_mul_ps(tW, m1->r[3]);
    tX = _mm_add_ps(tX, tZ);
    tY = _mm_add_ps(tY, tW);
    tX = _mm_add_ps(tX, tY);
    res.r[3] = tX;

    return res;
}
/* Quaternion Computation Operations */
static GH_INLINE GHvec ghQuatMul(GHvec q0, GHvec q1) {
    /*
    x = y0 * z1 - z0 * y1 + x0 * w1 + w0 * x1;
    y = z0 * x1 - x0 * z1 + y0 * w1 + w0 * y1;
    z = x0 * y1 - y0 * x1 + z0 * w1 + w0 * z1;
    w = w0 * w1 - x0 * x1 - y0 * y1 - z0 * z1;
    */
    GHvec t0, t1, r;
    t0 = _mm_shuffle_ps(q0, q0, _MM_SHUFFLE(3, 0, 2, 1));
    t1 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(3, 1, 0, 2));
    r  = _mm_mul_ps(t0, t1);
    t0 = _mm_shuffle_ps(q0, q0, _MM_SHUFFLE(0, 1, 0, 2));
    t1 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(0, 0, 2, 1));
    t0 = _mm_mul_ps(t0, t1);
    r  = _mm_sub_ps(r, t0);
    t0 = _mm_shuffle_ps(q0, q0, _MM_SHUFFLE(1, 2, 1, 0));
    t1 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(1, 3, 3, 3));
    t0 = _mm_mul_ps(t0, t1);
    t0 = _mm_mul_ps(t0, nGHNegateW.v);
    r  = _mm_add_ps(r, t0);
    t0 = _mm_shuffle_ps(q0, q0, _MM_SHUFFLE(2, 3, 3, 3));
    t1 = _mm_shuffle_ps(q1, q1, _MM_SHUFFLE(2, 2, 1, 0));
    t0 = _mm_mul_ps(t0, t1);
    t0 = _mm_mul_ps(t0, nGHNegateW.v);
    r  = _mm_add_ps(r, t0);
    return r;
}

#endif
